module finite_diff_mod
  
  use constants_mod
  use namelist_mod
  use boundary_mod, only : pbc_idx

  implicit none

  private
  public :: compute_interface_flux_finite_diff
  public :: compute_flux_jacobian_finite_diff
  
  contains

  subroutine compute_interface_flux_finite_diff(q, dt, flux)

    real(r8), intent(in) :: q(1-nhalo:nx+nhalo)
    real(r8), intent(in) :: dt
    real(r8), intent(out) :: flux(0:nx)
    integer i

    select case (trim(space_scheme))
    case ('finite_diff_2nd')
      do i = 0, nx
        flux(i) = u * 0.5_r8 * (q(i+1) + q(i))
      end do
    case ('finite_diff_4th')
      do i = 0, nx
        flux(i) = u / 12.0_r8 * (7.0_r8  * (q(i+1) + q(i)) - &
                                           (q(i+2) + q(i-1)))
      end do
    case default
      print*, 'Unknown space_scheme: ', trim(space_scheme)
      stop
    end select
    
  end subroutine compute_interface_flux_finite_diff

  subroutine compute_flux_jacobian_finite_diff(var, dt, jac)

    real(r8), intent(in) :: var(1-nhalo:nx+nhalo)
    real(r8), intent(in) :: dt
    real(r8), intent(out) :: jac(1:nx, 1:nx)
    integer i, idx_m1, idx_p1, idx_m2, idx_p2

    jac = 0.0_r8
    select case (trim(space_scheme))
    case ('finite_diff_2nd')
      do i = 1, nx
        idx_m1 = pbc_idx(i-1, nx)
        idx_p1 = pbc_idx(i+1, nx)
        jac(i, idx_m1) = -0.5_r8 * u
        jac(i, i     ) = 0.0_r8
        jac(i, idx_p1) = 0.5_r8 * u
      end do
    case ('finite_diff_4th')
      do i = 1, nx
        idx_m2 = pbc_idx(i-2, nx)
        idx_m1 = pbc_idx(i-1, nx)
        idx_p1 = pbc_idx(i+1, nx)
        idx_p2 = pbc_idx(i+2, nx)
        jac(i, idx_m2) = u / 12.0_r8
        jac(i, idx_m1) = -u / 12.0_r8 - u / 12.0_r8 * 7.0_r8
        jac(i, i     ) = 0.0_r8
        jac(i, idx_p1) = u / 12.0_r8 * 7.0_r8 + u / 12.0_r8
        jac(i, idx_p2) = -u / 12.0_r8
      end do
    case default
      print*, 'Unknown space_scheme: ', trim(space_scheme)
      stop
    end select
    jac = -jac / dx

  end subroutine compute_flux_jacobian_finite_diff

end module finite_diff_mod